﻿using DataAnalyticsTools.Models;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace DataAnalyticsTools.Core
{
    /// <summary>
    /// 趋势分析器
    /// 提供时间序列数据的趋势分析、变化点检测和季节性分解
    /// </summary>
    public static class SeasonalAnalyzer
    {

        /// <summary>
        /// 季节性趋势分解
        /// 将时间序列分解为季节性成分和趋势成分
        /// 使用加法模型：原始数据 = 趋势 + 季节性 + 残差
        /// </summary>
        /// <param name="data">单变量时间序列数据</param>
        /// <param name="period">季节性周期（必须 >= 2）</param>
        /// <param name="method">分解方法：MA(移动平均)或LOESS(局部回归)</param>
        /// <returns>包含趋势、季节性和残差成分的分析结果</returns>
        /// <example>
        /// 输入:
        ///   data = [2, 3, 4, 3, 2, 3, 4, 3]  // 周期为4的季节性数据
        ///   period = 4
        /// 输出:
        ///   SeasonalComponents = [-0.5, 0.5, 1.5, 0.5, -0.5, 0.5, 1.5, 0.5]  // 季节性模式（中心化）
        ///   TrendComponents = [2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5]        // 平稳趋势
        ///   Residual = [0, 0, 0, 0, 0, 0, 0, 0]                               // 残差
        /// </example>
        public static SeasonalAnalysisResult SeasonalTrendDecomposition(
            float[] data,
            int period,
            DecomposeMethod method = DecomposeMethod.MA)
        {
            // 输入验证
            if (data == null || data.Length == 0)
            {
                return new SeasonalAnalysisResult
                {
                    TrendComponents = new float[0],
                    SeasonalComponents = new float[0],
                    Residual = new float[0]
                };
            }

            if (period < 2)
                throw new ArgumentException("周期必须大于等于 2", nameof(period));

            if (period > data.Length)
                throw new ArgumentException($"周期 ({period}) 不能大于数据长度 ({data.Length})", nameof(period));

            if (data.Length < period * 2)
                throw new ArgumentException($"数据长度至少需要是周期的 2 倍（当前: {data.Length}, 需要: {period * 2}）", nameof(period));

            int n = data.Length;

            // 1. 提取趋势成分
            float[] trend = method == DecomposeMethod.MA
                ? ExtractTrendMA(data, period)
                : ExtractTrendLOESS(data, period);

            // 2. 去趋势化（移除趋势成分）
            float[] detrended = new float[n];
            for (int i = 0; i < n; i++)
            {
                detrended[i] = data[i] - trend[i];
            }

            // 3. 提取季节性成分
            float[] seasonal = ExtractSeasonal(detrended, period);

            // 4. 计算残差（原始数据 - 趋势 - 季节性）
            float[] residual = new float[n];
            for (int i = 0; i < n; i++)
            {
                residual[i] = data[i] - trend[i] - seasonal[i];
            }

            return new SeasonalAnalysisResult
            {
                TrendComponents = trend,
                SeasonalComponents = seasonal,
                Residual = residual
            };
        }

        #region 季节分解核心方法

        /// <summary>
        /// 使用移动平均法提取趋势
        /// 对于偶数周期，使用 2×period 的中心化移动平均
        /// 对于奇数周期，使用标准的中心化移动平均
        /// </summary>
        private static float[] ExtractTrendMA(float[] data, int period)
        {
            int n = data.Length;
            float[] trend = new float[n];

            // 对于偶数周期，使用经典的季节性分解方法
            // 需要先计算 period 的移动平均，再计算 2×1 的移动平均
            if (period % 2 == 0)
            {
                // 第一步：计算 period 的移动平均
                float[] ma1 = new float[n];
                int halfPeriod = period / 2;

                for (int i = 0; i < n; i++)
                {
                    int start = Math.Max(0, i - halfPeriod);
                    int end = Math.Min(n - 1, i + halfPeriod - 1);

                    float sum = 0;
                    int count = 0;
                    for (int j = start; j <= end; j++)
                    {
                        sum += data[j];
                        count++;
                    }
                    ma1[i] = count > 0 ? sum / count : data[i];
                }

                // 第二步：计算 2×1 的移动平均（中心化）
                for (int i = 0; i < n; i++)
                {
                    if (i == 0)
                    {
                        trend[i] = ma1[i];
                    }
                    else if (i == n - 1)
                    {
                        trend[i] = ma1[i];
                    }
                    else
                    {
                        trend[i] = (ma1[i - 1] + ma1[i]) / 2.0f;
                    }
                }
            }
            else
            {
                // 奇数周期：使用标准的中心化移动平均
                int halfPeriod = period / 2;

                for (int i = 0; i < n; i++)
                {
                    int start = Math.Max(0, i - halfPeriod);
                    int end = Math.Min(n - 1, i + halfPeriod);

                    float sum = 0;
                    int count = 0;
                    for (int j = start; j <= end; j++)
                    {
                        sum += data[j];
                        count++;
                    }

                    trend[i] = count > 0 ? sum / count : data[i];
                }
            }

            return trend;
        }

        /// <summary>
        /// 使用LOESS(局部加权回归)提取趋势
        /// LOESS = Locally Weighted Scatterplot Smoothing
        /// 在每个点周围使用加权最小二乘法进行局部线性回归
        /// </summary>
        private static float[] ExtractTrendLOESS(float[] data, int period)
        {
            int n = data.Length;
            float[] trend = new float[n];

            // 带宽参数：控制局部回归的窗口大小
            // 通常取周期的1.5-2倍，但至少要有足够的点进行回归
            float bandwidth = Math.Max(period * 1.5f, Math.Min(n * 0.3f, 15));

            // 最小带宽保护
            const float MIN_BANDWIDTH = 3.0f;
            bandwidth = Math.Max(bandwidth, MIN_BANDWIDTH);

            for (int i = 0; i < n; i++)
            {
                // 计算权重
                float[] weights = new float[n];
                float sumWeights = 0;

                for (int j = 0; j < n; j++)
                {
                    float distance = Math.Abs(i - j);

                    // 使用三次权重函数（Tricube kernel）
                    if (distance < bandwidth)
                    {
                        float u = distance / bandwidth;
                        // w(u) = (1 - |u|^3)^3 for |u| < 1
                        float cubicU = u * u * u;
                        float weight = (1 - cubicU);
                        weight = weight * weight * weight;
                        weights[j] = Math.Max(0, weight);
                        sumWeights += weights[j];
                    }
                    else
                    {
                        weights[j] = 0;
                    }
                }

                // 防止除零错误
                if (sumWeights < 1e-10f)
                {
                    trend[i] = data[i];
                    continue;
                }

                // 归一化权重
                for (int j = 0; j < n; j++)
                {
                    weights[j] /= sumWeights;
                }

                // 执行加权线性回归（局部线性拟合）
                // 计算加权均值
                float weightedMeanX = 0;
                float weightedMeanY = 0;
                for (int j = 0; j < n; j++)
                {
                    weightedMeanX += weights[j] * j;
                    weightedMeanY += weights[j] * data[j];
                }

                // 计算斜率
                float numerator = 0;
                float denominator = 0;
                for (int j = 0; j < n; j++)
                {
                    float dx = j - weightedMeanX;
                    float dy = data[j] - weightedMeanY;
                    numerator += weights[j] * dx * dy;
                    denominator += weights[j] * dx * dx;
                }

                // 计算趋势值
                if (Math.Abs(denominator) > 1e-10f)
                {
                    float slope = numerator / denominator;
                    float intercept = weightedMeanY - slope * weightedMeanX;
                    trend[i] = slope * i + intercept;
                }
                else
                {
                    // 如果无法计算斜率，使用加权平均
                    trend[i] = weightedMeanY;
                }
            }

            return trend;
        }

        /// <summary>
        /// 提取季节性成分
        /// 通过对每个季节位置的去趋势数据求平均来识别重复的季节性模式
        /// 最后进行中心化处理，使季节性成分的平均值为0
        /// </summary>
        /// <param name="detrended">去趋势后的数据</param>
        /// <param name="period">季节性周期</param>
        /// <returns>季节性成分数组（与输入等长，平均值为0）</returns>
        private static float[] ExtractSeasonal(float[] detrended, int period)
        {
            int n = detrended.Length;
            float[] seasonal = new float[n];

            // 第一步：计算每个季节位置的平均值
            // 例如：对于周期=4的数据，计算索引 0,4,8,... 的平均值作为第1个季节位置的值
            float[] seasonalPattern = new float[period];
            int[] counts = new int[period];

            for (int i = 0; i < n; i++)
            {
                int seasonIndex = i % period;
                seasonalPattern[seasonIndex] += detrended[i];
                counts[seasonIndex]++;
            }

            // 第二步：计算每个季节位置的平均值
            for (int i = 0; i < period; i++)
            {
                if (counts[i] > 0)
                {
                    seasonalPattern[i] /= counts[i];
                }
                else
                {
                    // 理论上不应该发生，但作为保护
                    seasonalPattern[i] = 0;
                }
            }

            // 第三步：中心化季节性成分（使平均值为0）
            // 这确保了：原始数据 = 趋势 + 季节性 + 残差
            // 其中趋势和季节性都不会相互影响
            float seasonalMean = seasonalPattern.Average();
            for (int i = 0; i < period; i++)
            {
                seasonalPattern[i] -= seasonalMean;
            }

            // 第四步：将周期性模式扩展到整个时间序列
            for (int i = 0; i < n; i++)
            {
                seasonal[i] = seasonalPattern[i % period];
            }

            return seasonal;
        }

        #endregion
    }
}
